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We describe an analysis of the First International Pulsar Timing Array Data Challenge, which 
was designed to test the ability of new and existing algorithms to constrain the properties of a 
stochastic gravitational-wave background influencing the arrival-times of pulsar signals. We employ 
a robust, unbiased Bayesian framework developed by van Haasteren to study the three Open and 
Closed datasets of the IPTA data-challenge. We test various models for each dataset and use 
MultiNest to recover the evidence for the purposes of Bayesian model-selection. The parameter 
constraints of the favoured model are confirmed using an adaptive MCMC technique. Our results 
for ClosedI favoured a gravitational-wave background with strain amplitude at / = 1 yr -1 , A, 
of (1.1 ± 0.1) x 10~ 14 , power spectral-index 7 = 4.30 ± 0.15 and no evidence for red-timing noise 
or single-sources. The evidence for Closed2 favours a gravitational-wave background with A — 
(6.1±0.3) x 10~ 14 , 7 = 4.34±0.09 with no red-timing noise or single-sources. Finally, the evidence 
for Closed3 favours the presence of red-timing noise and a gravitational-wave background, with no 
single-sources. The properties of the background were A = (5±1) x 10~ 15 and 7 = 4.23±0.35, while 
the properties of the red-noise were Af re< j = (12±4) ns and 7 rc d = 1.5±0.3. In all cases the redness of 
the recovered background is consistent with a source-population of inspiraling supermassive black- 
hole binaries. We also investigate the effect that down-sampling of the datasets has on parameter 
constraints and run-time. Finally we provide a proof-of-principle study of the ability of the Bayesian 
framework used in this paper to reconstruct the angular correlation of gravitational-wave background 
induced timing-residuals, comparing this to the Hellings and Downs curve. 

PACS numbers: 



I. INTRODUCTION 

There is a large international effort focussed towards 
the first direct detection of gravitational waves (GWs). 
The existing and planned ground-based instruments, 
such as AdLIGO [J, AdVirgo and KAGRA [3], are 
kilometre-scale interferometers sensitive to frequencies 
^10— 10 3 Hz, and will likely be operating at design sen- 
sitivity by the end of this decade. In addition, it is hoped 
that a space-based interferometer with arm-lengths of 
~ 10 9 m, such as eLISA/NGO 4J, will be operable by 
the end of 2020s, and sensitive in the range ~ 0.1 — 100 
mHz. 

The first indirect confirmation of the existence of GWs 
came from precision timing of the pulsar PSR B1913+16 
[5], whose inferred binary orbital-energy loss was found 
to be consistent with the prediction of general relativ- 
ity. This precision analysis was made possible by the of- 
ten sub-/zs level of timing precision achieved through the 
measurements of pulsar-signal time-of-arrivals (TOAs) 
[5J [7] whose accuracy can rival that of atomic clocks. 

The precision of millisecond pulsars can be exploited 
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for GW detection through the use of pulsar timing arrays 
(PTAs) [8 1, which can effectively use the Milky Way as 
a kpc-scale GW-detector. Tens of Galactic millisecond 
pulsars have been observed over several years to search 
for the influence of a GW perturbing the space-time met- 
ric along each pulsar-Earth line-of-sight . PTAs are 
sensitive to low GW frequencies (1 — 10 nHz), where this 
range is set by the observational time-span (/i ow ~ 
and the cadence (/high ~ 1/(2AT)), and as such are com- 
plementary to other GW detection experiments. 

It is not only GWs which can induce deviations of the 
TOAs from a timing model. The dominant perturbation 
is caused by the deterministic spin-down of the pulsar 
itself, as its rotational energy is extracted to power the 
EM outflow. There are also stochastic contributions to 
the deviations caused by a variety of sources, including 
clock noise, receiver noise and variations of the dispersion 
measure of the intervening interstellar medium. These ef- 
fects must be accounted for and removed from the TOAs 
to produce the timing residuals, which then contain only 
the influence of unmodelled phenomena, including GWs. 
While there is a rich literature on the subject of the de- 
tection of single GW-sources using pulsar-timing (e.g., 
[IME]), the most likely source of GWs for PTAs is a 
stochastic gravitational- wave background (GWB). 

An isotropic, stochastic GWB may be created by a 
superposition of many single sources which are not indi- 
vidually resolvable. In the PTA band, the largest contri- 
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bution will likely come from a cosmological population 
of inspiraling supermassive-black-hole-binary (SMBHB) 
systems, with typical masses ~ 10 4 — 10 10 M Q . 

The fractional energy-density of the Universe in a GW- 
background is usually given as, 



«GW(/) = 



1 dp GW {f) 
Pc d(\nf) 



(1) 



where / is the observed GW- frequency, p c = 3H 2 /8ttG 
is the energy-density required for a fiat Universe, and 
h c (f) is the characteristic strain of the GW-background 
in a frequency interval centred at /. 

The characteristic strain spectrum of a GW- 
background resulting from inspiraling binary systems 
is approximately h c (f) oc /~ 2 / 3 [19-22 . We can ap- 
proximate the characteristic strain spectrum of a GW- 
background from other sources as a power-law also. Some 
measurable primordial background contributions may 
have a power-law index of —1 [23l [24], while the back- 
ground from decaying cosmic strings [25H28| may have 
— 7/6 [20] . For most models of interest, we can describe 
an isotropic, stochastic GW-background by [50] . 



h c (f) = A 



J_ 
yr - 



(2) 



This characteristic strain spectrum is related to the one- 
sided power spectral density of the induced timing resid- 
uals by, 
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where 7 = 3 — 2a. 

Hcllings and Downs [31j developed a simple cross- 
correlation technique for pulsars affected by the same 
stochastic, isotropic GWB, showing that the cross- 
correlation of the induced timing-residuals has a distinc- 
tive angular signature dependant only on the angular- 
separation of the pulsars: 

3 , . , 1 11 
Cab = ^ ln ( x ) ~ 4* + 2 + 2 ' ^ ' 

where x = (1 — cos# a b)/2, and 6 ao is the angular separa- 
tion of the pulsar sky-locations. 

While a GWB induces correlated residuals which typi- 
cally have a steep, red spectrum, there may be additional 
uncorrelated red-noise contributions from rotational ir- 
regularities in each individual pulsar [32] , for which the 
spectrum is given as, 



S{f) = N 2 cd 
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Over the last several years constraints on the ampli- 
tude of an isotropic GWB have been published by the 
three major PTA collaborations [55H55] . the European 



Pulsar Timing Array (EPTA) [SB], the North Amer- 
ican Nanohertz Observatory for Gravitational Waves 
(NANOGrav) [57], and the Parkes Pulsar Timing Ar- 
ray (PPTA) [55]. The International Pulsar Timing Ar- 
ray (IPTA) [50] consortium combines these three efforts 
and recently initiated the first "IPTA Data Challenge" 
[40) whose aim was to test new and existing algorithms 
for the purpose of constraining the properties of a back- 
ground of GWs using PTAs. 1 

In this paper we describe an analysis of the first IPTA 
data challenge. We use the time-domain Bayesian frame- 
work developed by van Haasteren et al. [53J,[H]. Al- 
though faster implementations of this method have re- 
cently been suggested [42] [43], we employ the original 
framework and analyse the uncompressed data for the 
purposes of model-selection with the recovered Bayesian 
evidence. To evaluate the evidence of posterior param- 
eter distributions, we use the MultiNest nested sam- 
pling algorithm, and confirm favoured-model parameter 
constraints with an adaptive MCMC algorithm. 

This paper is organised as follows. In Sec. [Tl] we de- 
scribe Bayesian inference and how the evidence, which is 
the key quantity for model-selection, is computed. Sec- 
tion III describes pulsar timing analysis, including how 
the raw data is processed and a derivation of the likeli- 
hood of the timing-residuals given the timing-model and 
GWB parameters. The stochastic sampling techniques 



we have employed are described in Sec. IV followed by 
a description of the IPTA challenge data in Sec.jVlRe- 
sults are presented in Sec. [yTj followed in Sec. |VII| by 



brief studies of acceleration of the algorithm via down- 
sampling of the data and using the data to recover the 
Hellings and Downs correlation curve. We finish with our 
conclusions in Sec. IVIIII 



II. BAYESIAN INFERENCE 

Bayes' theorem states that the posterior probability 
density function (PDF), p(/Z|D, of the parameters /2 
describing a hypothesis model H, and given data D is 

p(D\p,M)p(ix\H) 



(6) 



P(D\H) ' 
where, 

p(D\fl,H) = C(fT) — likelihood of data given parameters, 
p{p\'H) = 7r(/7) = prior PDF of parameters, 
p(D\H) = Z = Bayesian evidence. (7) 

The Bayesian evidence, Z, is the probability of the 
observed data given the model % 



(8) 



1 Full details of this first challenge will be made available in a 
forthcoming IPTA data-challenge paper. 
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TABLE I: An interpretation of the Bayes factor in determining 
which model is favoured, as given by Jeffreys 1441 . 



Bayes factor, K. 


ln(/C) 


Strength of evidence 


< 1 : 1 


< 


Negative (supports H\) 


1 : 1 to 3 : 1 


0-1.1 


Barely worth mentioning 


3 : 1 to 10 : 1 


1.1 - 2.3 


Substantial 


10 : 1 to 30 : 1 


2.3 - 3.4 


Strong 


30 : 1 to 100 : 1 


3.4-4.6 


Very strong 


> 100 : 1 


> 4.6 


Decisive 



For posterior inference within a model, Z plays the role of 
a normalisation constant and can be ignored. However, 
if we want to perform model selection then this evidence 
value becomes key. In Bayesian model comparison we 
compute the Bayes factor 

p(n 2 \D) = p0\H 2 )p(H 2 ) = Z 2 x p(H 2 ) (Q) 
pCH^D) p0\Ui)p(Ui) -ZixpW 

where p^H^) / pi^-i) is the prior probability ratio for the 
two competing models. This can often be set to one, and 
we will do so in the remainder of this analysis. The Bayes 
factor is then just the evidence ratio. Since the evidence 
is the average of the likelihood over the prior volume, it 
automatically incorporates Occam's razor, which states 
that, all else being equal, a model with fewer parameters 
is favoured. Hypothesis Hi is chosen if the Bayes factor is 
sufficiently large. Jeffreys |Hj gave a scale interpretation 
for the Bayes factor, which is shown in Table [Tj 

We employ Bayesian model-selection in the following 
study to determine which phenomena provide the best 
explanation for the observed pulsar TOAs in the first 
IPTA data challenge. 

III. PULSAR TIMING ANALYSIS 

Observations of pulsars lead to measurements of the 
pulsar TOAs. The emission-time of a pulse is given in 
terms of the observed TOA by [13 EH] , 

C^^-Aq-Ats-Ab, (10) 

where A Q is the transformation from the site TOAs 
to the Solar-system barycentre, Ajg accounts for the 
delaying-effects as the pulse propagates through the in- 
terstellar medium, and Ab converts to the pulsar-frame 
for binary pulsars. 

In the first IPTA data challenge [40] the raw data is 
in the form of pulsar parameter files (".par") and timing 
files (".tim"). The parameter file contains first estimates 
of the pulsar timing-model parameters; these parameters 
describe deterministic contributions to the arrival times. 
The vector of measured arrival times will be composed 
of a deterministic and a stochastic contribution (from 



time-correlated stochastic signals which are modelled by 
a random Gaussian process), 

The stochastic process has auto-correlation, 

Ca = (5trstf p ), (12) 

where the elements of the covariance matrix are 
parametrised by a set of parameters, 4>. Using the 
Wiener-Khinchin theorem, we can then define the auto- 
correlation as the Fourier transform of the power spectral 
density, 

/>oo 

CK0= / S(f) C os(f nj )df, (13) 
Jo 

where r,j = 27r|tj — tj\, and S(f) is the power spectral 
density of the time-series Stf sp . A closed- form expres- 
sion for the auto-correlation of a time-series influenced by 
an underlying power-law PSD is given in van Haasteren 
et al. 35], and is used in the following. 

A. Processing raw arrival-times 

The ".par" and ".tim" files are fed to the Tempo2 
software package [45H47] which processes the raw arrival- 
times. A vector of "pre-fit" timing-residuals are com- 
puted using the first guesses, /3o,i, of the "to" timing- 
model parameters from the ".par" files. This first guess 
is usually precise enough so that a linear approximation 
can be used in the TOA fitting procedure, so that the 
post-fit timing residual are 

5t = 6&> l1 + Af£ (14) 

where SiP lf are the pre-fit timing-residuals (length n), 
£ is the vector of deviations from the pre-fit parameters 
(length to) defined as £ a = f3 a — /3 , a , and M is the (tixto) 
"design-matrix" , describing how the residuals depend on 
the timing-model parameters. Tempo2 does not take 
into account the possible time-correlated stochastic sig- 
nal in the TOAs, but performs a weighted least-squares 
fit for the timing-model parameter values. Hence it is 
possible that some of the time-correlated stochastic sig- 
nal is absorbed in this fitting procedure, which is unde- 
sirable. 

The Tempo2 analysis provides output-residuals and 
the design matrix, M. The design matrix describes the 
dependence of the timing residuals on the timing-model 
parameters. The output-residuals form the input data 
vector for further study. 

B. Generalised least-squares (GLS) estimator of 
stochastic and deterministic parameters 

We now want to use the Tempo2 output-residuals 
to determine any correlated stochastic signal affecting 
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the pulse arrival times. We assume that the part of 
the stochastic signal removed by the fitting procedure is 
small, so that the Tempo2 output-residuals are related 
linearly to the stochastic contribution to the residuals 



St = ^' sp + Af£ 



(15) 



where C = C^-C^M (A^C^M) M t C-\ When 
dealing with large datasets and many pulsars, C involves 
the multiplication and inversion of high dimensional ma- 
trices. In the case of multiple pulsars, the design matrix, 
covariance matrix and total residual vector are given by, 



where, in this case, St refers to the output-residuals from 
Tempo2. We note that the £ appearing in this equation 
is different from that appearing in Eq. ( 14 ) . 

The stochastic timing residuals, 5t ISP , arise from a 
time-correlated stochastic process with covariance ma- 
trix C. This covariance matrix may contain contribu- 
tions from the GWB, white-noise from TOA-errors, and 
possibly red-timing noise which is uncorrelated between 
different pulsars. The likelihood of measuring post-fit 
residuals, St, given the fit parameters £ and stochastic 
parameters, </>, is, 



£(sM, 



i 



^(27r)"detC ' 
1 



exp 



St 



Nltf C- 1 (St-M£) 



(16) 



This likelihood expression is effectively a GLS estimator, 
and is the basis for the framework used in this paper to 
study the first IPTA data challenge. 

If we assume flat priors on the timing-model 
parameters then these parameters can be analyti- 
cally marginalised over. The posterior distribution 
marginalised over timing-model parameters is |48) . 
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(18) 



where C a f, is the auto-covariance matrix between pul- 
sars a and b, M a are the individual pulsar design matri- 
ces and St a are the individual pulsar residual vectors. We 
can split C a b into contributions from various stochastic 
sources. So, the covariance between the ith residual of 
pulsar a and the jth residual of pulsar b is, 



C, 



(ai)(b 3 ) 



- °(oi)(6j)" 



^TOA 
-°(oi)(&2)- 



^EQUAD 
~°(oi)(«) " 



(19) 



(17) where, 
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oo 

r(i - 7rcd ) sm (^) (f inj r^ - £(-ir 



(2n)!(2n+l- 7rod ) 



(20) 



In the order listed, these are i) the GW-background 
covariance, CSuy) , which depends on £ a b, the Helling- 
Downs correlation between pulsars a and &; ii) the TOA 
error-bar covariance, C™^^, arising from white noise in 
each individual pulsar, which is characterised by a sep- 
arate, pre-specified and fixed, amplitude Ai( oi ) for each 



pulsar a and time i, plus an overall scaling factor, EFAC, 
which is common to all pulsars but is allowed to vary 
as a model parameter; iii) the covariance, C^^^ , of 
an additional white-noise which is common to all pulsars 
and time-stamps and characterised by a single amplitude 
parameter EQUAD; and iv) the covariance, C^ )(bj y of 



5 



a red timing-noise in each pulsar, which is modelled as 
a power-law with amplitude iV re( j and slope 7 rc d- The 
timing-model fit reduces the sensitivity of the residuals 
to the low- frequency cutoff, which is required when 
7 > 1. Hence, provided fiT << 1 we can ignore all 
terms with n > 2 in the infinite summation. To avoid 
numerical artefacts we choose fi = 10 -3 • r ' 



Expression (17) can be written more compactly and in 
a way which is slightly faster to compute [H]: 

P($\5t) = x 



exp 



^/(27r)"-™dct(G T C*G) 
-SPG (G T CGy 1 G T Si 



(21) 



where G is the matrix constructed from the final (n — m) 
columns of the matrix U in the SVD of the design matrix, 
AI = UYjV* . The matrix G can be pre-computed and 
stored in memory for use in each likelihood calculation. 



Equation (21) provides a robust, unbiased Bayesian 



framework for the search for correlated signals in PTAs, 
and is used in all of the following analysis. 



where A.hij = hij(t p ,Q) — hij(t e ,fl), is the difference in 
the metric perturbation at the pulsar and at the solar 
system barycentre. We ignore the uncorrelated pulsar- 
term in this analysis. 

This frequency-shift is integrated over time to give the 
induced timing residuals, 



R(t) 



z(t')dt'. 



(24) 



We will search for monochromatic and burst sources and 
descriptions of the models used to describe these sources 
are given in Appendix [A"] 



IV. STOCHASTIC SAMPLING TECHNIQUES 

We now discuss two different stochastic sampling tech- 
niques used to reconstruct the posterior PDF of the 
model parameters and to calculate the evidence value 
for Bayesian model-selection. 



C. Including single-sources in the search 

While this analysis focuses on the detection and char- 
acterisation of a background of GWs, we will also con- 
sider the possible presence of a single monochromatic 
or burst source perturbing the arrival times of pulses. 
In practice, it may be necessary to include single- 
source models in all background searches to allow the 
background-induced residuals to be described by Gaus- 
sian statistics [H5] . There is a large literature on the sub- 
ject of single-source detection in the context of PTAs. In 
particular, Sesana et al. |T7j describe how the commonly 
used approximation of a single power-law spectrum for a 
stochastic background from inspiraling SMBHBs breaks 
down at frequencies higher than 10 -8 Hz due to the dom- 
inance of single sources. Likewise, the authors of [13] - 
[TS1 [TB] have studied the ability of a PTA to infer the 
presence of multiple resolvable monochromatic sources, 
and to constrain their properties. 

We use the formalism of [48] to combine the search for 
single sources with the search for a background. The for- 
malism is a simple modification to the background search, 
in which the residuals are now described by, 

5t = (5?' sp + s + M£, (22) 

where s is the deterministic contribution to the residuals 
from a single source. 

Models for the +,x GW polarisation amplitudes can 
be used to compute the frequency-shift of pulses induced 
by the GW. The redshift of signals from a pulsar in the 
direction of unit vector p, induced by the passage of a 
GW coming from direction f2 is, 

^ 1 4- \l ■ n 



A. Markov chain Monte Carlo (MCMC) sampling 

Markov chain Monte Carlo (MCMC) techniques pro- 
vide an efficient way to explore a model parameter space. 
An initial point, x~o, is drawn from the prior distribution 
and then at each subsequent iteration, z, a new point, 
y, is drawn from a proposal distribution, q(y\x) and the 
Metropolis-Hastings ratio evaluated, 

A random sample is drawn from a uniform distribution, 
u G U[0, 1], and if u < R the move to the new point is 
accepted and we set afj+i = y . If u > R, the move is 
rejected and we set Xi+\ = xl. 

The MCMC samples can be used to carry out integrals 
over the posterior 

1 N 

i— 1 

The ID marginalised posterior probability distributions 
in individual model parameters then follow by binning 
the chain samples in that parameter. 

The trick to using this technique efficiently is to choose 
an appropriate proposal distribution. In our analysis we 
employ an adaptive MCMC procedure, which utilises an 
'in-flight' estimation of the sampled-chain's covariance 
matrix to construct an updating proposal distribution. 
This covariance matrix is updated at each iteration, with 
a certain chain memory [50H52] , We use several of the 
procedures outlined in 52 . A full description of this 
technique can be found in Taylor and Gair [531 and ref- 
erences therein]. 
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A single likelihood calculation on one core can take 
as much as ~ 40 seconds, so it is very important that 
we achieve a fast burn-in. We do this for each dataset 
by using the built-in threading capability of LAPACK 
[53] . This allows as many as 12 cores to speed-up a single 
matrix multiplication operation, and ultimately reducing 
the likelihood calculation time to ~ 5 seconds. Hence we 
use 5 independent chains with threading to burn-in. Af- 
ter approximately 2 hours burn-in is achieved, and we can 
sample from the end of these chains to initiate a larger 
run of 512 cores, with no threading, to collect posterior 
samples. Collecting ~ 10 5 samples takes approximately 
3 hours. 



B. Nested Sampling & MultiNest 

The nested sampling algorithm is a Monte Carlo 
method, originally proposed by Skilling [55] for evaluat- 
ing the Bayesian evidence, Z. For a full description of the 
MultiNest algorithm see Feroz and Hobson [55], Feroz 
et al. [5 7) . but we describe the basics in the following 
section. 

The basic idea is to populate parameter space with 
"live" points drawn from the prior. These points move as 
the algorithm proceeds, climbing together through nested 
contours of increasing likelihood. At each iteration, the 
points are ordered in terms of their likelihood, and the 
point with lowest likelihood is replaced by a point with 
higher likelihood than this lowest-likelihood point. 

The biggest difficulty in nested sampling is to effi- 
ciently sample points of higher likelihood to allow the 
live-points to climb. If we were to simply draw points 
from the prior volume, then the acceptance rate of new 
points in the live-set would steadily decrease, since at 
later iterations the live-set occupies a smaller and smaller 
volume of the prior space as it climbs. MultiNest over- 
comes this drawback by using a sophisticated ellipsoidal 
rejection-sampling technique, whereby the current live- 
set is enclosed by (possibly overlapping) ellipsoids, and 
a new point drawn uniformly from the enclosed region. 
This technique successfully copes with multimodal dis- 
tributions and parameter spaces with strong, curving de- 
generacies. 

The evidence is calculated by transforming the multi- 
dimensional integral in Eq. Q into a one-dimensional 
integral which is easily numerically evaluated. We define 
the prior volume, X as, 

dX = n(jl)d N (m, (27) 

such that, 

X(X) = f n(fl)d N fx, (28) 
Jc(p)>\ 

where the integral extends over the region of the N- 
dimensional parameter space contained within the iso- 



likelihood contour C{p) = A. Hence, Eq. ^ can be writ- 
ten as, 

Z= [ CdX, (29) 
Jo 

where C(X) is a monotonically decreasing function of X. 
If we order the X values (0 < X M < . . . < X 1 < X = 1), 
then the evidence, Z can be approximated numerically 
using the simple trapezium rule, 

M 

Z = ^C iWi , (30) 

i=l 

where the weights, Wi, are given by Wi = 
{Xi-i — Xi + i) /2. 

As a by-product of the exploration of the parameter 
space by the evolving live-set, MultiNest also permits 
reconstruction of the parameter posterior PDFs. Once Z 
is found, the final live-set, as well as the discarded points, 
are collected and assigned probability weights to give the 
posterior probability of each point. These points can be 
binned to give full and marginalised posterior PDFs. 

With MultiNest 's built-in MPI routines we use ~ 
800 live-points for all runs, and typically employ 160 
cores such that MultiNest finishes in less than 12 hours. 



V. DESCRIPTION OF DATA 

The first IPTA data challenge consists of three "Open" 
and three "Closed" datasets. For the "Open" datasets, 
the properties of the noise and background were given, 
allowing the calibration and testing of algorithms. The 
parameters of the "Closed" datasets are not due to be 
revealed until after the deadline. However the format of 
the data is the same in both sections. In all data sets, we 
have 36 pulsars distributed across the sky, each having 
130 pulse arrival times measured over an average time 
span of 5 years. 

A. Open data 

The open section of the challenge consisted of three 
separate sets of data, increasing in the level of complexity. 
All three data sets contained a GWB background with 
a power spectral density slope given by 7 = 3 — 2a = 
13/3, which is consistent with a background induced by 
a population of inspiraling SMBHBs. Each data set was 
for a total time span of 5 years with 130 observations per 
pulsar. The data sets differed in their sampling cadence 
and in the TOA noise in the pulsars. 

In "OpenI" the intrinsic noise in the pulsar TOAs 
was white, with the same amplitude of 100 ns in each 
pulsar, and the sampling cadence was uniform with a 
rate of one sample every two weeks. The characteristic 
strain-spectrum amplitude of the GWB at / = 1 yr _1 
was A = 5 x 10~ 14 . 
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In "Open2" the intrinsic TOA noise was again white, 
but the amplitude was different for each pulsar, with nom- 
inal values as given in Table [TTJ These white noise levels 
are consistent with realistic levels measured for IPTA 
pulsars. In addition, the sampling rate was no longer 
uniform but random, with an average cadence of 2 weeks 
±5 days. The characteristic strain-spectrum amplitude of 
the GWB was i = 5x 1CT 14 . 

In "Open3" the intrinsic pulsar noise has both a white 
component with levels as in "Open2" and an uncorre- 
cted red component, which had the same power spec- 
trum for each pulsar, but with a different realisation for 
each pulsar's TOAs. The red-noise power spectrum is 
S(f) = 5.77 x 10- 22 sec 1 ' 3 /- 1 - 7 , where / is in Hz. This 
is equivalent to the expression, 

«(fl = ^(j^i)( 1 ^i)" 7 " 4 , m 

where 7 rc d = 1.7 and N vc d = 10.1 ns. The characteristic 
strain-spectrum amplitude of the GWB was A = 10 -14 . 



VI. RESULTS 
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FIG. 1: The ID and 2D marginalised posterior distributions of 
OpenI. The star and dashed lines show the injected values of 
parameters. 



(half-width of 68% credible region) for each parameter. 
All plots were made with GetDist in the CosmoMC pack- 
age [55] , 



We used MultiNest to compute evidence values for 
model selection and confirmed parameter constraints for 
the favoured model with adaptive MCMC. Except in 
certain circumstances which we will discuss, we made 
no attempt to recover the individual white-noise lev- 
els of each pulsar, nor have we permitted the pulsars 
in Open/Closed3 to have individual red-noise ampli- 
tude/index parameters. Rather, where necessary, pul- 
sars share a global TOA error-bar scaling (GEFAC), ex- 
tra white-noise contribution (GEQUAD), red-noise am- 
plitude and red-noise index. However, all white-noise and 
red-noise is uncorrelated between different pulsars. 

From experience with the open datasets, we found that 
uniform priors on all GWB and red-noise parameters was 
reasonable. A test-run on OpenI with uniform priors on 
the log of the GWB-amplitude yielded similar results. 
The posterior distribution was narrow compared to the 
prior which covered eight decades in amplitude, which 
explains why the results of uniforms priors on ampli- 
tudes and log-amplitudes were similar. When fitting for 
the EQUAD or white-noise parameters, our prior was 
iogio (ctequad/wn) € [-11.0,-3.0]. To allow for the 
influence of prior choice on the evidence calculation, we 
only consider a model disfavoured if the evidence is lower 
than another model by at least Aln(Z) = 3, which was 
the typical difference seen when different prior choices 
were tested. However, if two models have similar evi- 
dence but one has many more unconstrained parameters, 
we naturally favour the simpler model. 

Tablcs lTTTl and lTvl show all tested models for all datasets 
with the associated evidence. The strength with which 
certain models are favoured can be gauged using Table 
[TJ We quote maximum-a-posteriori values and la errors 



A. Open Datasets 

1. OPENi 

This dataset provided the first test of the analysis 
pipeline. The time-stamp of each data sample and the 
white noise level in the pulsar were read in as input, along 
with the post-fit residuals. Although these were constant 
for all pulsars in OpenI, this implementation allowed 
easy generalisation to the other data sets. The timing- 
models of all pulsars converged after a single Tempo2 
fitting procedure. 

From Table IIHI it is clear that the evidence deci- 
sively favours the presence of a GWB. The properties 
of this background are A = (4.8 ± 0.2) x 10~ 14 and 7 
1.1 '2 ± 0.08. Figure [l] shows the ID and 2D marginalised 
posterior distributions for these parameters. These re- 
covered parameters are consistent with the injected val- 
ues. 



2. Open^ 

The differences between this dataset and OpenI were 
the random observation-cadence, and different white in- 
trinsic TOA errors in each pulsar. These provided no 
difficulties for the pipeline. 



From Table III it is clear that the evidence decisively 
favours the presence of a GWB, although the recovered 
background-amplitude showed a 2a deviation from the 
injected value. Adding in GEFAC and GEQUAD pa- 



TABLE II: The white noise levels of each pulsar in the 0pen2 and Open3 datasets [58] . 



Pulsar 


RMS WN (jjs) 


Pulsar 


RMS WN ( M s) 


Pulsar 


RMS WN ( M s) 


Pulsar 


RMS WN (/xs) 


J0030+0451 


0.31 


J1022+1001 


0.37 


J 1730-2304 


0.83 


J1910+1256 


0.17 


J0218+4232 


4.81 


J1024-0719 


0.25 


J1732-5049 


1.74 


J1918-0642 


0.87 


J0437-4715 


0.03 


J1045-4509 


2.68 


J1738+0333 


0.24 


J1939+2134 


0.02 


70^1 ^ 0900 


A ^ 


71/1^ 3 3^0 


1 £0 


717/11^-1^^1 
J 1 ( *±i-\- loOl 




Ti q^_l-9Q0& 


O 1 8 
U.lo 


J0621+1002 


9.58 


J 1600-3053 


0.23 


J1744-1134 


0.14 


J2019+2425 


0.66 


J071 1-6830 


1.32 


J 1603-7202 


0.70 


J1751-2857 


0.90 


J2124-3358 


1.52 


J0751+1807 


0.78 


J1640+2224 


0.19 


J1853+1303 


0.17 


J2129-5721 


0.87 


J0900-3144 


1.55 


J1643-1224 


0.53 


J1857+0943 


0.25 


J2145-0750 


0.40 


J1012+5307 


0.32 


J1713+0747 


0.04 


J 1909-3744 


0.04 


J2317+1439 


0.25 



TABLE III: Models tested for all Open datasets. The most 
favoured model for each dataset is shown in bold. The acronyms 
correspond to: Gravitational- Wave Background (GWB), 
Time-Of- Arrival (TOA), White Noise (WN), Global EFAC 
(GEFAC), Global EQUAD (GEQUAD) and Red-Noise (RN). 



Model 


ln(JZ) 


OpenI 


Null (TOA-crrors only) 


-190502 


GWB, TOA-errors 


61405.2 ± 0.1 


Open2 


Null (TOA-errors only) 


-821338 


GWB, TOA-crrors 


55466.2 ± 0.1 


GWB, GEFAC, GEQUAD, TOA-crrors 


55471.7 ± 0.1 


GWB, WN(J1857+0943) 


55574.8 ± 0.1 


GWB, WN(J0437-4715, J1857+0943, 


55568.2 ± 0.2 


J 1909-3744) 




WN(J1857+0943) 


-819198 


Open3 


Null (TOA-crrors only) 


37379 


GWB, TOA-crrors 


56038.9 ± 0.1 


GWB, RN, TOA-errors 


56037.8 ± 0.1 



rameters improved the evidence and reported GEFAC = 
1.06 ± 0.01 and GEQUAD consistent with zero. This 
was revisited after the analysis of Closed2 described 
below. Visual inspection of the output-residuals showed 
"glitchy" behaviour in J1857+0943, and reprocessing 
the arrival-times by Tempo2 indicated that the timing- 
solution had not converged after one fitting-iteration. 
Rather than re-fit we allowed the total white-noise in this 
pulsar to be a model parameter to be fitted. This created 
a huge improvement in the evidence and the recovery of 
the injected GWB parameters. 

Fitting for the white-noise of several other pulsars 
which showed "glitchy" behaviour shifted the max-a- 
posteriori background-amplitude slightly closer to the in- 
jected value, but the evidence was slightly worse. We 
thus conclude that J1857+0943 was the main source 
of the inconsistency between the injected and recovered 
background-amplitude. This was the only pulsar whose 
timing-solution did not converge after a single Tempo2 



fitting procedure. The recovered properties of the GWB 
are A = (5.4 ± 0.3) x 10" 14 and 7 = 4.33 ± 0.09, with 
an effective EFAC on J1857+0943 of 2.2 ± 0.2. Figure [2] 
shows the ID and 2D marginalised posterior distributions 
for these parameters. The recovered GWB properties are 
consistent with the injected values. 



3. Open 3 



This dataset had all of the complexity of Open2 with 
the addition of red timing-noise. This red-noise had a 
common amplitude and index for all pulsars, but was un- 
corrclatcd between different pulsars. The timing-models 
of all pulsars converged after a single Tempo2 fitting 
procedure. 



From Tabic III it is clear that the evidence decisively 
favours the presence of a GWB. The low- frequency be- 
haviour in the residuals was dominated by the GWB, 
such that in the model with a GWB and red timing-noise 
the red-noise properties remained largely unconstrained. 
The model with only a GWB gave slightly better evi- 
dence and the GWB parameter constraints were consis- 
tent with the GWB+RN model. We present results for 
the GWB+RN model since we know that that this is the 
true description of the dataset, but in a blind analysis we 
would not have inferred the presence of red-noise. The 
properties of the GWB are A = (1.17±0.14) x 10~ 14 and 
7 = 4.1 ± 0.2, which are consistent with injected values. 
Figure [3] shows the ID and 2D marginalised posterior dis- 
tributions for all parameters. The null model, where only 
TOA-errors are used to fit the timing-residuals, provides 
much greater evidence in this case than in OpenI /Open2 
since the injected GWB amplitude is five times smaller 
than in the previous two datasets. 
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FIG. 2: The ID and 2D marginalised posterior distributions of Open2. The star and dashed lines show the injected values of 
parameters. 



B. Closed Datasets 

1. ClosedJ 

A list of the models tested and the associated In Z 
values are in Table |IV| The timing-models of all pulsars 
converged after a single Tempo2 fitting procedure. We 
analysed ClosedI using the same method as OpenI. In- 
cluding a GWB and TOA errors provides evidence which, 
when compared to the null evidence, decisively proves 
the presence of a GWB. Comparing the evidence for a 
GWB to that for an uncorrelated red-noise process in 
each pulsar favours the correlated GWB. Red timing- 
noise in addition to the GWB is not strongly disfavoured, 
but did not improve the fit either. Adding in GEFAC 
and GEQUAD parameters did not improve the fit, as 
GEFAC was consistent with 1 and GEQUAD was con- 
sistent with zero. The model with a GWB, GEFAC and 
TOA-errors implied a GEFAC of 0.98 ± 0.01 i.e., a 2a 



deviation. However, this model did not significantly im- 
prove the fit. Including a monochromatic source did not 
improve the fit and all the parameters of the single source 
were unconstrained. 

We conclude that Closed 1 contains only a GWB with 
no red-noise or single GW sources. The properties of this 
GWB are A = (1.1 ± 0.1) x 10~ 14 and 7 = 4.30 ± 0.15. 
Figure [4] shows the ID and 2D marginalised posterior 
distributions for these parameters. 



2. Closed^ 

This dataset provided an interesting challenge. A list 
of the models tested and the associated InZ values are 
in Table El 

We found that assuming a model composed of a GWB 
and TOA-errors only did not provide a satisfactory fit. 
Including red timing-noise provided much higher evi- 
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FIG. 3: The ID and 2D marginalised posterior distributions of Open3. The stars and dashed lines show the injected values of 
parameters. 



dence, however the index of the red-noise reached the 
edge of it's prior at 7 rcc j = 1 (where the closed-form 
expression for the auto-covariance of a power-law PSD 
breaks down). A model with 7 rc( j = would effec- 
tively be a white noise model, which can also be in- 
cluded by the EQUAD parameter introduced earlier. The 
model with a GWB, GEFAC and GEQUAD provided 
the best evidence of these initial three models. This 
model implied A = (5.9 ± 0.5) x 10~ 14 , 7 = 4.45 ± 0.15, 
GEFAC = 0.89 ± 0.01 and GEQUAD = 289 ± 6 ns. 

We then tested for the presence of a single-source 
in the dataset. The evidence did appear to favour a 
single-source, either a monochromatic source or a GW 
burst. However, the recovered parameters favoured a 
source with a frequency several orders of magnitude 
above Nyquist, and a very well constrained sky-location 
close to the pulsar J0437-4715. As a sanity check, we 
removed this pulsar from the dataset and again tested 



for the presence of the source. In this case, all source 
properties were unconstrained, which is inconsistent with 
the presence of a single source, since the 35 remaining 
widely separated pulsars should provide adequate trian- 
gulation. Furthermore, the GEFAC and GEQUAD pa- 
rameters were then consistent with 1 and respectively. 
Pulsar J0437-4715 has arrival-times with nominally mea- 
sured timing precision of 0.03 /is. A visual inspection of 
the post-fit residuals showed that such small error-bars 
were not sufficient to explain the high-frequency fluctua- 
tions that were present in this pulsar in addition to the 
red-noise induced by the GWB. We reprocessed the ar- 
rival times through several Tempo2 iterations, which in- 
dicated that the residuals had not converged after a single 
fitting iteration. 

We therefore allowed the total white-noise in J0437- 
4715 to be a parameter to fit in our analysis. This cre- 
ated a very significant improvement in the evidence. Fit- 
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FIG. 4: The ID and 2D marginalised posterior distributions of 
ClosedI. 



TABLE IV: Models tested for all Closed datasets. The most 
favoured model for each dataset is shown in bold. The acronyms 
correspond to: Gravitational- Wave Background (GWB), 
Time-Of- Arrival (TOA), White Noise (WN), Global EFAC 
(GEFAC), Global EQUAD (GEQUAD) and Red-Noise (RN). 



Model 


ln(.Z) 


ClosedI 


Null (TOA-errors only) 


51584.5 


GWB, TOA-errors 


62028.1 ± 0.1 


RN, TOA-errors 


62015.8 ± 0.1 


GWB, RN, TOA-errors 


62026.7 ± 0.1 


GWB, RN, GEFAC, GEQUAD, TOA-errors 


62022.7 ± 0.1 


GWB, GEFAC, TOA-errors 


62025.2 ± 0.1 


GWB, Monochromatic-source, TOA-errors 


62028.3 ± 0.1 


Closed2 


Null (TOA-errors only) 


-815557 


GWB, TOA-errors 


53807.4 ± 0.1 


GWB, RN, TOA-errors 


54149.9 ± 0.1 


GWB, GEFAC, GEQUAD, TOA-errors 


54178.3 ± 0.1 


GEFAC, GEQUAD, TOA-errors 


50082.4 ± 0.1 


GWB, GEFAC, GEQUAD, Burst-source, 


54235.0 ± 0.2 


TOA-errors 




GWB, GEFAC, GEQUAD, Monochromatic-source, 


54222.8 ± 0.2 


TOA-errors 




GWB, WN(J0437-4715) 


55246.3 ± 0.1 


GWB, WN(J0437-4715, J1455-3330, 


55235.2 ± 0.2 


J1741 + 1351, J1909-3744) 




GWB, WN(J0437-4715), Monochromatic-source, 


55245.9 ± 0.1 


TOA-errors 




GWB, WN(J0437-4715), RN, TOA-errors 


55245.4 ±0.1 


RN, WN(J0437-4715) 


55227.3 ± 0.1 


WN(J0437-4715) 


-240920 


CLOSED3 


Null (TOA-errors only) 


53886.6 


RN, TOA-errors 


56168.9 ± 0.1 


GWB, RN, TOA-errors 


56185.6 ± 0.1 


GWB, RN, GEFAC, GEQUAD, TOA-errors 


56180.3 ± 0.1 


GWB, RN, Monochromatic-source, TOA-errors 


56185.6 ± 0.1 



ting for the white-noise of several other pulsars which 
showed "glitchy" behaviour did not improve the fit, so 
we conclude that the anomalous effects were dominated 
by J0437-4715. This was the only pulsar whose timing- 
solution did not converge after a single Tempo2 fitting 
procedure. The insufficient white-noise in this pulsar was 
mimicking a single, high-frequency GW-source in the pul- 
sar's vicinity. Testing a model which fits for a GWB, 
the white-noise in J0437-4715 and a single-source left the 
source-parameters unconstrained, although the evidence 
was only slightly lower. 

Our final results for Closed2 were A — (6.1 ± 0.3) x 
1(T 14 , 7 = 4.34 ± 0.09. The effective EFAC on J0437- 
4715 was 30 ± 2, but this was entirely due to the poor 
timing-model fit in this pulsar. Figure [5] shows the ID 
and 2D marginalised posterior distributions for these 
parameters. When we performed a final analysis of 
this data, where the residuals had been cycled through 
Tempo2 multiple times to ensure convergence of the 
timing-model fits, we found GWB parameter constraints 
consistent with these results, and again no compelling 
evidence for a monochromatic source. 



3. Closed,? 

For the analysis of Closed3 we used the same model 
as Open3, including a red timing-noise component. A 
list of the models tested and the associated InZ values 



are in Table IV The timing-models of all pulsars con- 
verged after a single Tempo2 fitting procedure. 

The evidence decisively favours the presence 
of red timing-noise and a GWB. A model with 
GEFAC+GEQUAD gave values consistent with 1 and 
0, respectively. A model with a monochromatic-source 
left the source-parameters unconstrained. We conclude 
that Closed3 contains a GWB with red-noise and 
no single GW sources. The properties of the GWB 
were A = (5 ± 1) x 10~ 15 and 7 = 4.23 ± 0.35, while 
the properties of the red-noise were A rc d = (12 ± 4) 
ns and 7 rcc i = 1.5 ± 0.3. Figure [6] shows the ID 
and 2D marginalised posterior distributions for these 
parameters. 



VII. ACCELERATION BY DOWN-SAMPLING 

The algorithm used in this work is computationally 
expensive and two recent proposals have been made 
of approaches to PTA data analysis that are faster, 
van Haasteren [32] used a high-fidelity data-compression 
technique, which utilises an interpolation scheme of the 
compressed covariance matrix elements to constrain A 
and 7. The Fisher information is used to determine how 
much of the data might be redundant. Lentati et al. 
[43] proposed a method that avoids the dense-matrix 
multiplications and inversions altogether by rephrasing 
the likelihood in terms of matrix-vector operations and 
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TABLE V: OpenI parameter Icr errors and run-times on a single Sandy-Bridge node with 16 cores, for varying fractions of the data. 



Fraction of data Cadence / weeks Parameter constraints Approx run-time 







A 1 xl0~ 14 


7 




0.1 


20 


2.04 


0.48 


1 min 


0.125 


16 


0.65 


0.21 


2.5 min 


0.167 


12 


0.37 


0.15 


5 min 


0.2 


10 


0.32 


0.13 


8 min 


0.333 


6 


0.25 


0.11 


37 min 


0.5 


4 


0.22 


0.097 


2 hr 


1.0 


2 


0.18 


0.084 


11.5 hr 



banded-matrix inversions. This is achieved by mod- 
elling the effect of the GWB on the timing residuals di- 
rectly in the time-domain by a small number of inde- 
pendent Fourier components, effectively eliminating the 
off-diagonal components with i ^ j in the GW correla- 
tion matrix, and avoiding an a-priori prescription for the 



GWB spectrum. 

Using a reduced number of frequency components is 
analogous to down-sampling the data in the time do- 
main, so we investigated a very simple down-sampling 
of the data in OpenI to determine if reasonable param- 
eter constraints could be derived when using a subset 
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FIG. 6: The ID and 2D marginalised posterior distributions of CL0SED3. 
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FIG. 7: OpenI parameter la errors and run-times on a single Sandy-Bridge node with 16 cores, for varying fractions of the data. 
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FIG. 8: The 2D posterior distributions with 68% and 95% confidence intervals for batches of data with 10%, 20% and 50% of the total 
OpenI dataset. Stars show the injected parameters. 



TABLE VI: OpenI MultiNest 


run-times for the full dataset 


with varying numbers of computer 


cores. 


Cores 


MultiNest run-time 


16 


11 hr 30 min 


160 


2 hr 30 min 


768 


hr 45 min 



of the data. We would expect that a GWB induces a 
low-frequency variation of the timing-residuals, such that 
short cadence observations are redundant. 

Tableland Fig.[7]show the parameter constraints ob- 
tained using down-sampled datasets comprising different 
fractions of the total data. The downsampling is accom- 
plished by increasing the length of time between observa- 
tions, by selecting every n'th time sample. The timing- 
model is first derived using the entire dataset then the 
subset of timing residuals to be used in the analysis are 
selected. The appropriate G matrix is derived by keep- 
ing only the corresponding rows of the design matrix and 
deleting all other rows. 

We see that, at least for the case of OpenI, using 
only 20% of the full OpenI dataset (corresponding to 
a cadence of 10 weeks) is sufficient to achieve parame- 
ter constraints which are comparable to constraints from 
the full dataset, in approximately ~ 1% of the time. The 
posterior widths AA and A7 are respectively ~ 80% and 
~ 50% wider than those obtained from an analysis of the 
full dataset. If each 20% subset of data contained the 
same information as the full dataset, we would expect 
the relative decrease in performance going from 100% to 
20% of the data to be ~ y/E, which would mean a ~ 125% 
increase in posterior width. We do better than this ex- 
pectation since we are discarding high frequency infor- 
mation, while the gravitational wave background power 
is predominantly at lower frequencies. We do not do 
much better than this, however, although the analysis is 
much quicker. The corresponding 2D posterior distribu- 
tions for batches of data with 10%, 20% and 50% of the 
total OpenI dataset are shown in Fig. [8j 

For completeness, in Fig. [9] we show the result of an 



analysis of a 20% batch of the Open3 data, comparing 
the GWB and intrinsic red-noise parameter constraints 
to the analysis of the complete dataset. Once again we 
see that, as a first-cut analysis of the data, this down- 
sampling technique allows us to narrow our search space 
when following-up with a full-dataset analysis, while also 
achieving comparable parameter constraints. However, 
much more work is required to determine the optimal 
way to down-sample in the more realistic case where we 
may have long gaps in data-taking. 

We also show in Table |VI| the speed-up achieved when 
more cores are used to exploit the built-in MPI routines 
in MultiNest. This speed-up is clearly not linear, so 
the appropriate compromise between computational ex- 
penditure and run-time is left at the discretion of the 
user. 



A. Stacking the down-sampled posteriors 

Ideally we would like to use all of the data in our anal- 
ysis, requiring a method to stack the posterior distribu- 
tions from each batch of data. Unfortunately, each batch 
of data is not an independent sample since the post-fit 
residuals are the result of a Tempo2 fit to the entire 
dataset. Furthermore, even without pre-fitting a timing- 
model to the entire dataset, if we were to combine pos- 
teriors from analyses of separate batches of raw-arrival 
times we would still effectively be assuming that each 
batch contained GWB-induced residuals drawn from a 
different realisation of the GWB-spectrum, which is not 
the case. 

Modulo these issues, it is still worthwhile to investi- 
gate the posterior widths and parameter bias that are 
obtained when using the naive approach of taking the 
posterior distributions obtained from analysing (n — 1) 
distinct subsets comprising 1 / n'th of the data as a prior 
for an analysis of the remaining n'th batch of data. We 
compared the result of using the posterior distributions 
from the analyses of 4 batches of 20% of the OpenI data 
as priors on the final batch to the full coherent analysis of 
the OpenI dataset, and this is shown in Fig.[l0] In this 
case, the stacking of the down-sampled posteriors pro- 
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(a) (b) 

FIG. 9: The 2D posterior distributions with 68% and 95% confidence intervals for a 20% batch of Open3 data (red) is compared to the 
results of a full-dataset analysis (blue). The left panel shows GWB parameters, and the right panel shows the intrinsic red- noise 
parameters. Stars show the injected parameters in each case. 
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FIG. 10: The 2D posterior distribution from a full coherent 
analysis of the OpenI dataset (blue) is compared to using the 
posterior distributions from the analyses of 4 batches of 20% of 
the OpenI dataset as priors on the analysis of the final batch 
(red). Credible regions are 68% and 95%, and the injected 
parameters are indicated by the star. 
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FIG. 11: The reconstructed angular-correlation of the isotropic, 
stochastic GWB in 20% of the OpenI data. We parametrise the 
correlation at scp = {1°, 30°, 60°, 100°, 140°, 180°}. A 
cubic-spline is then used to interpolate the correlation at all 
pulsar angular separations between 1° and 180°, with the 
correlation at 8 Bcp = 0° fixed at unity. The Hcllings and Downs 
curve is shown as a black line, the max-a-posteriori correlation 
function is shown as a dotted line and the envelope of splines in 
the 95% credible region of this reconstruction is shown in red. 



duces a narrower final posterior, since by assuming each 
batch of data is an independent sample we gain a factor 
of \f% by combining them, while as pointed out above 
the individual 20% posteriors are less than a factor of 
yfh wider than the posterior for the full dataset. The 
difference is not that large and the bias in the maximum 
a-posteriori values appears to be quite small. Most im- 
portantly, the stacked posterior remains consistent with 
the injected parameters in this case. For these reasons, 
and the fact that the analysis took ~ 5% of the run-time 
of the full analysis, we advocate this technique as a useful 
first-cut analysis of PTA data. More studies are required 
to understand under what circumstances the stacked pos- 
terior will no longer be consistent with the injected val- 
ues. 



B. Recovering the angular-correlation function 

The acceleration of the likelihood-evaluation due to 
down-sampling allows us to perform a novel test. If we 
have an isotropic, stochastic GWB with general relativity 
as the correct description of gravity then the angular cor- 
relation of GWB-induced residuals will be the Hellings 
and Downs curve. However if we have more than the 
usual two general relativity polarisations then the angu- 
lar correlation could deviate from this curve [60l IS] . We 
use 20% of the OpenI dataset to reconstruct the angular 
correlation of the GWB-induced timing-residuals. 

We replace the angular-correlation function ^ a ;, by, 
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FIG. 12: The reconstructed angular-correlation of the isotropic, 
stochastic GWB in 20% of the OpenI data. We parametrise the 
angular correlation with the same functional form as the Helling 
and Downs curve, but vary the numerical co-efficients as model 
parameters. The Hellings and Downs curve is shown as a black 
line, the max-a-posteriori correlation function is shown as a 
dotted line and the envelope of correlation functions in the 95% 
credible region of this reconstruction is shown in red. 
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FIG. 13: The recovered spectrum of the GWB from the 
analysis of 20% batches of the OpenI dataset. We parameterise 
S(f) at fT = {0.5, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0}, where T = 5 yr. We 
compute the autocovariance in the time-series by performing a 
simple trapezoidal integration over S(f) cos(/r;j), where 
r = 27r|tj — tj\. As expected, we have very little sensitivity to 
frequencies below 1/T, however the recovered spectrum is 
consistent with the injected power-law spectrum. 
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(32) 



where ci, C2, C3, C4, C5, cq are allowed to vary as model pa- 
rameters, with prior range € [—0.5,1.0]. A cubic-spline 
interpolation of the correlations between 1° and 180° is 
then used to calculate the correlation at all other pulsar 
angular-separations. The result of this analysis is shown 



in Fig. 11 where the reconstructed angular correlation 
is shown to be consistent with the Hellings and Downs 
curve. 

An alternative technique to probe the angular correla- 
tion is to assume the correlation has the same functional 
form as the Hellings and Downs curve, but with different 
numerical coefficients, 



Cab = P0% ln(x) + piX +p 2 + (1 - p 2 )5 a b, 



(33) 



where x = (1 — cos# scp )/2 and po;Pi;P2 are varied as 
model parameters, with prior ranges chosen to be sym- 
metric around the true Hellings and Downs values. The 
numerical coefficients for the Hellings and Downs curve 
are po = 1.5, p\ = —0.25 and p 2 — 0.5, so we chose pa- 
rameter prior ranges of po € [0.0,3.0], p\ € [—0.5,0.0], 
Pi G [0.0,1.0] respectively. However, we verified that 
our results were not dependent on the exact width of 



these priors by checking the posterior distributions was 
not influenced by the prior restrictions. The result of 
this analysis is shown in Fig. |12| where the reconstructed 
angular correlation is again shown to be consistent with 
the Hellings and Downs curve. 



Using the ansatzes in Eq. ( 32 ) and Eq. ( 33 ) we have re- 



constructed the angular-correlation of the GWB without 
massively expanding the dimensionality of our parameter 
space, and the fact that we have demanded smooth vari- 
ation of the angular-correlation reduces the susceptibility 
of our constraints to intrinsic pulsar noise-processes. The 
technique which exploits the ansatz in Eq. ( 32 1 is more 



general, where we have performed a model-independent 
reconstruction of the angular-correlation. 

This is a proof-of-principle that illustrates how a 
Bayesian framework can be used to infer that any GWB 
present in PTA data is correlated with the distinctive 
angular signature expected in general relativity. Fur- 
ther work is required to quantify how the precision of 
the correlation reconstruction depends on the quality of 
the data, the number of pulsars etc. and to explore its 
sensitivity to non-GR polarisation states and anisotropy 
in the background. 



C. Model-independent probes of the GWB 
spectrum 

Finally, rather than assume that the GWB spectrum 
is described by a power-law, we could parameterise the 
power-spectral density, £(/), at certain frequencies to 
probe the spectrum in a model-independent way [43] . 
One would expect this approach to fail and simply re- 
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turn the characteristic oc f~ 2 spectrum associated with 
spectral leakage that results from windowing a finite 
time-series, since the underlying background spectrum 
is steeper than f~ 2 . However, this model-independent 
technique was shown to be successful in 43J and we be- 
lieve that spectral-leakage is avoided since the model fit- 
ting by Tempo2 removes the low-frequency power in the 
time-series which would dominate the leakage signature. 
Meanwhile, the information about the background is pre- 
served through our knowledge of the design-matrix, M, 
which encodes information about the fitting procedure. 

Using batches of 20% of the OpenI dataset, we re- 
peated the stochastic background search but replaced 
the closed-form expression for the autocovariance of a 
time-series induced by a power-law background by a sim- 
ple trapezoidal integration over S(f) cos(frij) (exploiting 
the Wiener-Khinchin theorem), where r = 27r|tj — tj\ and 
S(f) is parametrised at certain frequencies. Since we are 
investigating a steep, red spectrum we parameterise S(f) 
at fT = {0.5, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0}, where T is the ob- 
servation span. The results are shown in Fig. 13 , where 
we see that, as expected, we have very little sensitivity 
to frequencies below 1/T, and the recovered spectrum is 
consistent with the injected power-law spectrum. 



VIII. CONCLUSIONS 

We have used a Bayesian time-domain method to com- 
pute solutions to the first International Pulsar Timing 
Array data challenge, providing constraints on the prop- 
erties of the injected isotropic, stochastic gravitational- 
wave background. The posterior probability distribu- 
tion of the stochastic background parameters is analyt- 
ically marginalised over all deterministic pulsar timing- 
model parameters. We tested this algorithm on the Open 
datasets, successfully recovering the values of the injected 
parameters. 

We made use of the MultiNest algorithm to calcu- 
late the Bayesian evidence value for various models and 
to provide posterior parameter distributions. Posterior 
PDFs were also obtained using an adaptive MCMC code 
to provide a cross-check of the results. Various mod- 
els were tested for each dataset, including null models 
where only TOA error-bars were present. The computed 
evidence values were then used for model-selection, de- 
termining which collection of astrophysical sources pro- 
vided the best explanation for the observed pulsar TOA 
deviations. 

The results for the Closed datasets were as follows. 
The evidence for ClosedI favoured a gravitational-wave 
background with strain amplitude at / = 1 yr _1 , A, of 
(1.1 ± 0.1) x 10~ 14 , spectral-index 7 = 4.30 ± 0.15 and 
no compelling evidence for uncorrelated red timing-noise 
or single-sources. The evidence for Closed2 favoured a 
gravitational-wave background with A = (6.1 ± 0.3) x 
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7 



4.34 ± 0.09, and, again, no compelling ev- 



idence for red timing-noise or single-sources. Finally, 



the evidence for Closed3 favoured the presence of red 
timing-noise and a gravitational-wave background, with 
no single-sources. The properties of the background were 
A = (5 ± 1) x KT 15 and 7 = 4.23 ± 0.35, while the 
properties of the red-noise were N re ^ = (12 ± 4) ns and 
7red = 1-5 ± 0.3. The adaptive MCMC sampling proce- 
dure provided consistent results in all cases. 

We found it was necessary to cycle the output-residuals 
through the Tempo2 software package (which performs a 
weighted least-squares fit of a deterministic timing-model 
to the pulse arrival-times) several times for some pulsars. 
Timing-model fits should be inspected to ensure conver- 
gence has been attained, otherwise a poor fit in even a 
single pulsar can severely bias any parameter constraints. 

Following our analysis of the full challenge datasets, we 
investigated the effect that down-sampling of the datasets 
has on parameter constraints and an the analysis run- 
time. This was achieved by only selecting every nth 
residual in the data set and design matrix. Given that the 
bottle-neck steps of the likelihood evaluation are 0(n 3 ) 
matrix-matrix operations, down-sampling provided a sig- 
nificant speed-up in the algorithm. When tested on the 
OpenI data, we found that using 20% of the full data 
could achieve comparable parameter constraints to the 
full dataset, but in only ~ 1% of the time. Stacking the 
posterior distributions of (n — 1) batches of data as pri- 
ors for the analysis of the remaining nth batch of data 
resulted in a small bias in the posterior width, but this 
was not severe and the posterior remained consistent with 
the injected parameters. 

As a final test, we used 20% of the OpenI data to in- 
vestigate whether we could constrain the angular corre- 
lation of the GWB-induced timing-residuals. If we have 
an isotropic, stochastic GWB, and general relativity is 
the correct description of gravity, then this angular cor- 
relation should be described by the Hellings and Downs 
curve, Eq. Q. We parametrised the correlation at vari- 
ous pulsar angular-separations, employing a cubic spline 
interpolation to compute the correlations at all other sep- 
arations. Doing so we found that the reconstructed angu- 
lar correlation was not consistent with zero or full corre- 
lation, but rather showed a distinctive angular signature 
which was consistent with the Hellings and Downs curve. 
This was the first test of a novel procedure, and signif- 
icantly more work is required to fully understand what 
can be accomplished in practice. We intend to further 
extend this line of inquiry, given that the determination 
of the angular correlation of GWB-induced residuals is 
key to understanding the effects of finiteness in the back- 
ground [49] and possible modifications to general relativ- 
ity [501 EE]. 

In our future work we will explore the optimal formal- 
ism to allow both background and single-source extrac- 
tion, as well as the effect that finiteness of the background 
has on inferred background parameter constraints and 
the signature of the angular correlation. This first IPTA 
data challenge will be followed up by more complex and 
realistic challenges designed to push current techniques 
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to their limit and encourage more optimal methods to 
be developed. We plan to analyse these future IPTA 
mock datasets, and real IPTA data, employing similar 
techniques as used here. 

Inferring the presence and properties of both 
a gravitational-wave background and individual 
gravitational-wave sources is a key aim of pulsar 
timing array projects. Data challenges such as this first 
IPTA challenge allows established members of the field 
and newcomers to test a variety of different approaches 
so that optimal techniques can be determined. By 
the end of the 2020's we are likely to have precision 
instrumental coverage over a wide range of GW frequen- 
cies, in the form of possible follow-ups to the advanced 
ground-based interferometers, the commissioning of a 
space-based interferometer, and the full operation of the 
SKA [62], which will be the most sophisticated radio 
telescope ever built. These complementary instruments 
will all contribute to the goal of gravitational-waves 
becoming a precision astronomical tool. 
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Appendix A: Single GW sources 
1. Monochromatic sources 

A monochromatic source in the PTA-band is likely to 
be an SMBHB in the very early inspiral stage. Sesana 
and Vecchio [TS] have shown that it is reasonable to ig- 
nore eccentricity, spin-effects and frequency-evolution for 
these systems. Hence the quadrupolar approximation for 
GW-emission can be used to describe the waveform. 

In our analysis, we use the residual signal model of 
Sesana and Vecchio [15]. The polarisation amplitudes 
for a monochromatic source of frequency / are, 

h + = A (1 + cos 2 i) cos ($(t) + $ ) , 

h x = „4costsin($(i) + $ ) . (Al) 

Hence for pulsar a, 

s a (t) = K [(1 + cos 2 i) F" (sin ($(*) + $ ) - sin $ ) 



where TZ = A/(2%f), and, 



F" = F C Q cos (2V>) + Ff sin (2tp) , 

F° = — F" sin (2-0) + Ff cos (2^) , (A3) 



where, 



F: = <j ^(sin 2 ( XQ )-2cos 2 ( Xtt ))sin 2 # 

+ 2 cos (x«) sin(x Q ) sin(20) cos(0 - j a ) 

+ cos 2 0) sin 2 ( XQ ) cos (20 - 2 la ) ' 



Ff = {- cos(x Q ) sin(xa) sin(0) sin(0 - j a ) 



1 



sin 2 (xa) cos(0) sin(27 Q - 20) 



l + h a ■ k 



1 



1 + n a ■ k 
(A4) 



and {7q, Xq} and {</>, 0} are the sky-locations of the pul- 
sar and GW-source in spherical polar co-ordinates, re- 
spectively. The azimuthal angle in spherical polars is 
equivalent to right ascension, while the polar angle is 
related to declination by = (ir/2 — dec). Further de- 
tails of the polarisation basis formalism can be found in 

G51I53I. 



2. Burst sources 

Various types of source could generate short-lived GW 
"bursts" in the PTA band, including cosmic strings, 
where cusps (sections of the string travelling at close to 
the speed of light) emit highly beamed radiation. Wave- 
form models of these cusp bursts exist, and take into 
account the spectrum of radiation when viewed slightly 
off the emission axis [Ml and references therein]. Bursts 
of gravitational radiation may also be emitted when two 
SMBHs pass close to one another on a highly eccentric or- 
bit, which can occur shortly after a major galactic merger 
[65] , However, for the purposes of identifying the pres- 
ence of a burst we adopt a simple sine-Gaussian waveform 
model, which is sufficiently generic to provide a good fit 
to many different burst sources. The waveform is cen- 
tred on a particular frequency, /o, and exponentially sup- 
pressed at nearby frequencies. In the time-domain this 
has the form l66ll67l. 



(2W (i-t b )) 2 



h+(t) = h +fi sin (27t/ < + $o) exp 
h x (*) = h xfi cos (2nf t + $ ) exp 



2Q 2 
(2nf (t-t b )f 



2Q 2 



(A5) 



+2 cos i (cos ($(£) + $ ) - cos $ 



(A2) 



As described before, the induced residuals at the Earth 
are given by modulating the metric perturbation by the 
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PTA antenna response function for each polarisation, and 
then integrating over time. A sine-Gaussian integrated 
over time is qualitatively similar to a sine-Gaussian. We 
therefore adopt the approach of _63J by searching for a 
residual-signal of sine-Gaussian form. Hence, 



s a (t) = R + F^ cos (2irf t + $ ) exp 
+ i? x F^sin(2^/ t + $ )exp 



(2nf Q (t-t b )Y 
2Q 2 

(27T/ (t - t h )f 



2Q 2 



(A6) 



where R+ = 1Z (l + cos 2 ij /2 and R x = 1Z cos / 



[1] G. M. Harry and the LIGO Scientific Collaboration, Clas- [28 
sical and Quantum Gravity 27, 084006 (2010). 
Virgo Technical Report VIR-0027A-09, Virgo (2009). [29 
K. Somiya and for the KAGRA Collaboration, ArXiv e- 

prints (2011), 1111 7185. [30 

eLISA/NGO, URL |http : //www . elisa-ngo . orgT] 
J. H. Taylor and J. M. Weisberg, ApJ 345, 434 (1989). 
V. M. Kaspi, J. H. Taylor, and M. F. Ryba, ApJ 428, 
713 (1994). [31 

D. N. Matsakis, J. H. Taylor, and T. M. Eubanks, A&A [32 
326, 924 (1997). 

R. S. Foster and D. C. Backer, ApJ 361, 300 (1990). [33 
M. V. Sazhin, Soviet Ast. 22, 36 (1978). 
S. Detweiler, ApJ 234, 1100 (1979). 

F. B. Estabrook and H. D. Wahlquist, General Relativity 
and Gravitation 6, 439 (1975). [34 
W. L. Burke, ApJ 196, 329 (1975). 

A. Sesana, A. Vecchio, and M. Volonteri, MNRAS 394, 
2255 (2009), 0809.3412. 

S. Babak and A. Sesana, Phys. Rev. D 85, 044034 (2012), [35 
1112.1075. 

A. Sesana and A. Vecchio, Classical and Quantum Grav- 
ity 27, 084016 (2010), 1001.3161. 

K. J. Lee, N. Wex, M. Kramer, B. W. Stappers, C. G. [36 
Bassa, G. H. Janssen, R. Karuppusamy, and R. Smits, | 
MNRAS 414, 3251 (2011), 1103.0115. [37 
A. Sesana, A. Vecchio, and C. N. Colacino, MNRAS 390, 
192 (2008), 0804.4476. [38; 
A. Petiteau, S. Babak, A. Sesana, and M. de Araujo, 
ArXiv e-prints (2012), 1210.2396. [39 
M. C. Begelman, R. D. Blandford, and M. J. Rees, Na- 
ture 287, 307 (1980). [40 

E. S. Phinney, ArXiv Astrophysics e-prints (2001), 
arXiv:astro-ph/0108028. [41 
A. H. Jaffe and D. C. Backer, ApJ 583, 616 (2003), 
arXiv:astro-ph/0210148. [42 
J. S. B. Wyithe and A. Loeb, ApJ 590, 691 (2003), [43 
arXiv:astro-ph/0211556. 

L. P. Grishchuk, Pis ma Zhurnal Eksperimental noi i Teo- [44 
reticheskoi Fiziki 23, 326 (1976). 

L. P. Grishchuk, Physics Uspekhi 48, 1235 (2005), 
arXiv:gr-qc/0504018. [45 
A. Vilenkin, Phys. Rev. D 24, 2082 (1981). 

A. Vilenkin, Physics Letters B 107, 47 (1981). [46 
S. Qlmez, V. Mandic, and X. Siemens, Phys. Rev. D 81, 
104028 (2010), 1004.0890. [47 



S. A. Sanidas, R. A. Battye, and B. W. Stappers, 

Phys. Rev. D 85, 122003 (2012), 1201.2419. 

T. Damour and A. Vilenkin, Phys. Rev. D 71, 063510 

(2005) , arXiv:hep-th/0410222. 

F. A. Jenet, G. B. Hobbs, W. van Straten, R. N. Manch- 
ester, M. Bailes, J. P. W. Verbiest, R. T. Edwards, A. W. 
Hotan, J. M. Sarkissian, and S. M. Ord, ApJ 653, 1571 

(2006) , arXiv:astro-ph/0609013. 

R. W. Hellings and G. S. Downs, ApJ 265, L39 (1983). 
R. M. Shannon and J. M. Cordes, ApJ 725, 1607 (2010), 
1010.4794. 

P. B. Demorest, R. D. Ferdman, M. E. Gonzalez, D. Nice, 
S. Ransom, I. H. Stairs, Z. Arzoumanian, A. Brazier, 
S. Burke- Spolaor, S. J. Chamberlin, et al., ArXiv e-prints 
(2012), 1201.6641. 

D. R. B. Yardley, W. A. Coles, G. B. Hobbs, J. P. W. 
Verbiest, R. N. Manchester, W. van Straten, F. A. Jenet, 
M. Bailes, N. D. R. Bhat, S. Burke-Spolaor, et al., MN- 
RAS 414, 1777 (2011), 1102.2230. 

R. van Haasteren, Y. Levin, G. H. Janssen, K. Lazaridis, 
M. Kramer, B. W. Stappers, G. Desvignes, M. B. Purver, 
A. G. Lyne, R. D. Ferdman, et al., MNRAS 414, 3117 
(2011), 1103.0576. 



European Pulsar Timing Array, URL http://www.epta. 
^-org/| 



North American Nanohertz Observatory for Gravita- 
tional Waves, URL http://nanograv.org/ 



Parkes Pulsar Timing Array, URL http : / /www . atnf . 



csiro . au/research/pulsar/ppta/ 

International Pulsar Timing Array, URL |http : / /www . | 
ipta4gw.org/ 



First IPTA Data Challenge, URL http : //www. ipta4gw. 
org/?page_id=214 

R. van Haasteren and Y. Levin, ArXiv e-prints (2012), 
1202.5932. 

R. van Haasteren, ArXiv e-prints (2012), 1210.0584. 

L. Lentati, P. Alexander, M. P. Hobson, S. Taylor, and 

S. T. Balan, ArXiv e-prints (2012), 1210.3578. 

H. Jeffreys, Theory of probability, International series of 

monographs on physics (Clarendon Press, 1983), ISBN 

9780198531937. 

G. B. Hobbs, R. T. Edwards, and R. N. Manchester, 

MNRAS 369, 655 (2006), arXiv:astro-ph/0603381. 

R. T. Edwards, G. B. Hobbs, and R. N. Manchester, 

MNRAS 372, 1549 (2006), arXiv:astro-ph/0607664. 

G. Hobbs, F. Jenet, K. J. Lee, J. P. W. Verbiest, D. Yard- 



20 



ley, R. Manchester, A. Lommen, W. Coles, R. Edwards, 
and C. Shettigara, MNRAS 394, 1945 (2009), 0901.0592. 
R. van Haasteren, Ph.D. thesis, Ph. D. thesis, University 
of Leiden (2011). pp. 176 (2011). 

V. Ravi, J. S. B. Wyithe, G. Hobbs, R. M. Shannon, 
R. N. Manchester, D. R. B. Yardley, and M. J. Keith, 
ArXiv e-prints (2012), 1210.3854. 

H. Haario, E. Saksman, and J. Tamminen, Compu- 
tational Statistics 14, 375 (1999), URL |http : //www . | 
springerlink.com/index/10. 1007/s001800050022 
H. Haario, E. Saksman, and J. Tamminen, Bernoulli 
7, pp. 223 (2001), ISSN 13 507265, URL |http : //www . | 
j stor . org/stable/3318737 

J. Dunkley, M. Bucher, P. G. Ferreira, K. Moodley, 
and C. Skordis, MNRAS 356, 925 (2005), arXiv:astro- 
ph/0405462. 

S. R. Taylor and J. R. Gair, Phys. Rev. D 86, 023502 

(2012), 1204.6739. 

LAPACK, URL |http : //www . netlib . org/lapack7| 
S. J., in Bayesian Inference and Maximum Entropy Meth- 
ods in Science and Engineering, edited by T. U. V. Fis- 
cher R., Preuss R. (Am. Inst. Phys., New York, 2004), 
vol. 735, p. 395. 

F. Feroz and M. P. Hobson, MNRAS 384, 449 (2008), 
0704.3704. 

F. Feroz, M. P. Hobson, and M. Bridges, MNRAS 398, 
1601 (2009), 0809.3437. 



URL http : //www. ipta4gw. 



URL http: //www. 



[58] Pulsar white-noise levels, 

| org/?page_id=148[ 

[59] CosmoMC, URL http://cosmologist.info/cosmomc/ 
[60] K. J. Lee, in American Institute of Physics Conference 
Series, edited by M. Burgay, N. D'Amico, P. Esposito, 
A. Pellizzoni, and A. Possenti (2011), vol. 1357 of Amer- 
ican Institute of Physics Conference Series, pp. 73-76, 
1105.5562. 

[61] S. J. Chamberlin and X. Siemens, Phys. Rev. D 85, 

082001 (2012), 1111.5661. 
[62] Square Kilometre Array, 

skatelescope . org/ 
[63] M. Pitkin, ArXiv e-prints (2012), 1201.3573. 
[64] F. Feroz, J. R. Gair, P. Graff, M. P. Hobson, and 

A. Lasenby, Classical and Quantum Gravity 27, 075010 
(2010), 0911.0288. 

[65] L. S. Finn and A. N. Lommen, ApJ 718, 1400 (2010), 
1004.3499. 

[66] B. Abbott, R. Abbott, R. Adhikari, J. Agresti, P. Ajith, 

B. Allen, R. Amin, S. B. Anderson, W. G. Ander- 
son, M. Arain, et al., Phys. Rev. D 77, 062004 (2008), 
0709.0766. 

[67] J. Abadie, B. P. Abbott, R. Abbott, T. D. Abbott, 
M. Abernathy, T. Accadia, F. Acernese, C. Adams, 
R. Adhikari, C. Affeldt, et al., Phys. Rev. D 85, 122007 
(2012), 1202.2788. 



